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ABSTRACT 


A 

A least squares estimator 4r is derived for the state transition 
matrix 4 of a linear, stationary sampled data system operating In a 

A 

stochastic environment. The estimator 4 Is shown to be unbiased and 
minimum varlanca under the condition of full observability of the state 
vector of the system. The estimator Is also shown to be the Maximum 
Likelihood Estimator for the case of the stochastic environment having 
Gaussian statistics. The estimation scheme Is compared with two 
other recently published estimation schemes, both of which are shown 
to be special cases of the scheme herein presented. 
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Vectors, (All lower case letters are vectors unless otherwise 
stated.) 

continuous representation of state vector 
discrete representation of state vector 
state vector of an equivalent system 
observation of the system 
system excitation 
additive measurement noise 

vector formed from a time sequence of a diuCrete state 
mode 

vector formed from a time sequence of an excitation mode 
i^ row vector in the matrix ® 


2. Matrices. (All upper case letters are matrices unless otherwise 

stated.) 

F matrix describing the homogeneous system 

D input matrix of the system 

4 (t) homogeneous solution of the system 

4 state transition matrix of the discrete representation 

r input distribution matrix of the discrete system 

J matrix containing the cost function 

X matrix formed from a collection of state vectors 

A 

4 estimate of the state transition matrix 

R finite approximation to the autocorrelation function 

P inverse of R 

m o 

H observability transformation 

4 * state transition matrix of an equivalent linear system 

Q covariance of the randcxn excitation 

3. Scalar quantities. 

L quadratic cost function 

2 

L((0i .0^) log of the likelihood function of a random sample 





n system order 

m sample size 

k time index 

a a scalar function defined in the recursive form of the 

identification algorithm - 


4. Operators. 

T 

X 

denotes x transpose 

(. 

denotes ( . ) inverse 

V-^ 

denotes gradient of ( . ) with respect to x 

tr( . ) 

denotes the trace of ( . ) 

E( . ) 

denotes the expectation or mean of ( . ) 

cov( . ) 

denotes the covariance of ( . ) 

(a,,) 

denotes A 
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1. • lie- ■ 
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‘I'qinji ' .. .,^1 w s- ^ I 
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1. INTRODUCTION 


The problem to be considered is the identification, in a stochastic 
environment, of linear, stationary (constant coefficient), sampled data 
process dynamics. The problem is more concisely defined as follows: 

Given: 

1) A system whose behavior may be described by a set of linear, 
constant coefficient, differential equations. 

2) A statistical description of the excitation of the system (the 
vector w in the figure below). 

3) A sequence of observations on the state vector (x^ in the 
figure). 

Problem: 

Determine a "best" estimate of the system dynamics. 



SYSTEM WITH 

i 

x(t) j 

w 

UNKNOWN 


DYNAMICS 



A. 


Figure 1. Schematic representation of the problem 

Identification, defined in this manner, amounts to estimating the 
solution of the differential equation postulated in 1) above. 

This problem may appear at the outset to be a rather restricted 
one. It becomes apparent, however, that a study of the problem takes 
on added significance when viewed in the context of a more general 
question. Given a vector of random time functions produced by some 
unkrxjwn process, can the process be adequately represented by a 
linear, stationary model, and if so, which linear model Is the best rep¬ 
resentation of the system? This work is motivated primarily with the 
hope of finding application in the investigatlpn of such questions. 

Significant contributions in this area have been made by Ho and 
^alen [ 2 3 and Lee L 3 ] . The purpose of this work will be to extend 
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the methods suggested by the above authors. The approach here will be 

• l i ( I 

to derive a scheme which minimizes a quadratic cost criteria, subject 
to some rather limiting assumptions. An attempt will be made to analyze 
the consequence of relaxing these restrictive assumptions. Finally, 
the methods of Ho and Whalen [ 2 ] and Lee [ 3 ] will be generated 
from the method herein derived. Acknowledgement is given to LT Ralph 
E. Hudson, USN, who first suggested the algorithm to be presented. 

2. BASIC MATHEMATICAL MODEL 

I 

Consider a dynamic system whose behavior may be defined by a 
set of n linear, constant coefficient, differential equations. 

X c F X + Dw* 

The solution is given by 

t 

x(t) ■ 4 (t-t ) X + J* 4(t-T) Dw*(t) d T 
0 0 * 

o 

I - 

where 4 (t-t ) satisfies 
o 

i 4 (t-g - F « (t-y 

1 . 

Introducing a sampling device of period T and a zero order hold 
on the excitation signal w*(t) makes possible the representation of the 
system at multiples of the sampling instant T by ' 

Vl ■ **k*rw^* 

• > .T.I ■ ■ . ; ! . I . ■.-■N 'i 

where 

£ Xj^ ■ x(kT) 

-r-i- N ''5 •• 'i. • 

4 ■ 4CF) 

. :iv:' ,.'j.: .:nh' T 

. u* * ■ rCD-J- ,4,(T-T)DdT, 

o 

u»!j 3 'u ■' '1 ii‘i ' . i . , 

p 03 -M: ' ‘ ‘ !'■ ■. Jh - ■•i i ■ j - . 
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Figure 2. Schematic representation of the difference equation 


3. DERIVATION OF THE IDENTIHCATION ALGORITHM 
Consider a process of the form 


Vi 


4 X, + w, 


(3.1) 


Here is defined to be a gausslan sequence of zero mean, with co- 
variance matrix Q. Suppose that a sequence of m+1 observations are 
made on the state vector x; (x^, x^, . . . x^, from which it is desired 
to form an estimate of the state transition matrix 4, assumed to be un- 

A 

known. The so-called "lease squares" estimator is that estimator 4 
which minimizes the quadratic cost function L (a scalar), where 


„ 'Vi - * *k>\+i - * 

k»“0 


(3.2) 


Equivalently, 4 minimises the trace of J, where 
n-1 o i A T 

<Vi - * V<*k+1 -« 

k«0 


(3.3) 


From Inspection, L * tr(J). 

ailL. I ■ , 

The classical method of solving a problem of this nature is to set 

A 

the gradient of the cost function with respect to the estimator 4 equal 

i ! 





to zero. The solution to this gradient equation Is sufficient to es¬ 
tablish a minimum cost. The method cannot be applied here, however. 

A 

since the gradient operation is not defined for the matrix . 

To find a solution, it is proposed to reformulate the minimization 

\ 

criteria. and then show the subsequent solution also minimizes the 
original cost function. tr(J). 

Let the order of the system be n. and partition the state transition 

T 

matrix 4 into n row vectors . 1*1. 2. . . . n. 

Then the process model 

Vi ■ 

may be written as m scalar equations. 

*k+i ‘ 

' .1 ' 

Here the superscript denotes the i^^ scalar component (or mode) 

of the vector. From the m+1 vector observations, (x . x,. . , . x ). 

0 1 m 

define 


CVi- • ' 

(nxm matrix) 

(mxl vector) (mxl vector) s. ... m 



From (3.5). these newly defined terms are related by 




(3.6) 


A ' ' n - T . 

Now consider an estimator , which minimizes the quadratic coet 

C 'J i MV'C , ■ J - . 

function g.(a scalar), where 

.1* _o ;rf' J T 'A(. ' < .. 


J U .T 




(3.7) 
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Comparing equations 3.3 and 3.7« It is seen that 


®1 " ^1 ' ^ * ^Ik^ 


Thus. the sum of the new cost function g^, i^l. 2. . . . n, Is pre¬ 
cisely the same as that previously defined. In that 
n 

I g. - tra) (3.8) 

^ 1-1 


A sufficient condition for the minimisation of tr(J). then, is that g. be 

A ^ 

minimised for 1-1, 2, . . . n. Further, g is not a function of 0. , 

y ^ A ^ 

l^J. This implies that the optimum estimator ^ does not depend upon 

. f ' ' * 

the choice of the estimator selected for the Jth row of the matrix es- 

A A 

tlmator, 4 . Thus, each 0^ may be found independently. 

To accomplish the minimisation, expand equation 3.7 and perform 
the gradient operation. 

’ . -vJx + 2 X^X (3.8) 

A 

Setting g^ - 0 and solving for 


- Vj’^X(X^X)"^ 


(3.9) 


T 

Equation (3.9) assumes X X to be nonsingular. 


The conditions for 


nonsingularity are discussed In Section 5. 

A 

The optimum estimator 4 can now be formed by placing the row 

A J 

vectors 0 ^ in matrix form. < 


•A 

« 



X(X^X)”^ 

n' 


(3.10) 











Finally, writing 4 as a function of the original sequence of vector 

observations x , x,, . . . x , gives 
o 1 m 


^ m-1 T r t”I 

*' ^ *k+l *k ^ *k*k - 

k*0 ^ ^ ^ k“0 ^ ^ 

Here £ x. x. is recognized as a finite approximation to the 
k-0 

autocorrelation function, R(kT), evaluated at k^O. Similarly, 


L ^ finite approximation to R(T). Making this identification 

k«0 


' ' i - R, (R (3.12) 

1 o 

■ » , 

where 

R, » R(T) R^ - R(0) . 

1 o 

4. DEVELOPMENT OF A RECURSIVE FORM 

A 

Note that (3.12) gives the optimum estimator 4 , based on a 
sample of size m+1. The estimator is, then, a function of m, allowing 

• A 

the functional notation 

! 0 ' ci A :• c ■ R, (m) [R (m)!”^ » - (4.1) 

m 1 o 

To keep the index consistent with the sample size, define ^ 

' A:" , 

m-1 ^ 

R,(m) ■ Z, X x^ r 

l Kfl K 

m-l 

R (m) ■ L X 

o k-0 

Incrementing the index in (4.1) gives 

' *■' ' r/ry ! 

^m+l ' tRj(m+l)][R^(m+l)r^ j J | 

A i 




From the definition of (m) 




To ease the reader through the remainder of the recursive development« 
let 

X Px X "x ,, R » R (m) 

m 1 m+1 o o 


Then 


Cl ■ (Ri+XixV„ + xxV‘ 


Appealing to the matrix Inversion lemma, [ 5 ] 


(R +xx^)“^ » r“^(I - xCx^ r“^ x+ 1]”^ x"^ r“^) 
o o o o 


(4.3) 


Recognizing R.r"^ as ^ 

1 o m 


'‘o ^f'"*** * ''o ^ 


‘ -1 -1 

Define the scalar a “ (x R^ + 1) 

m mom 


Collecting terms and simplifying 


m+1 m m m+1 m m m o 


(4.4) 


Define? »=R“^(m). 
m o 

From equations (4.3) and (4.4), the final recursive forms can be 
written down. 




^ +a (x .,-4 x)x^P 

m+1 m m m+1 m m m m 

' ^ ' 'Illy A ":; ' 


(4.5) 


u .;r:Jr.vg. y 'S'---;- 'u n ii? --r 


rir 
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(4.6) 


m+1 


P (I - tt X X P ) 
m m m m m 


a =(x P x^ + 1)^ 
m m m m 


(4.7) 


Implicit in the use of these recursive forms is the utilization of 

(4.1) to generate an initial estimator 4 and an initial matrix P , This 

m 

fact would seriously complicate a mechanization of the Iterative process. 

However, this difficulty can be avoided by the appropriate assignment 

of values to the initial estimator 4 , and the initial P matrix. Digi- 

m 

tal simulation has produced no apparent degradation of the estimation 
scheme when using this ruse. This subject is further discussed in 
Section 10. 

5. NECESSARY AND SUFFICIENT CONDITIONS FOR IDENTIFIABILITY 
Equation (3.9) assumed the nonsingularity of 
m-1 „ 


< «= E X, x: 

0 , . Tc Tc 

k*=0 


(nxn) matrix 


Consider the decomposition of *= X X 


o 

2 


n 


I X 


. X 


m-l 


n 

'm-l 


(nxm) matrix 


Here the superscript again denotes the element or mode of thq state 

vector. ' 

From the fact that X must be of rank n if R is to be non-singular, 

o 

it is immediately apiparent that m ^n is a necessary condition for 


1*4 




identlfiablllty. The number of observations on the state vector 

must be equal to or greater than the system order, n. From Inspection 

of X, it is seen that if any mode of the system, , k = 0, 1, . . . m-1, 

remains at zero for all k, or if any two modes remain at some constant 

value, X will be of rank r <n, and R will be singular. These con> 

° -1 

dltlons imply that all modes of the system must be excited if is 
to exist. 

Finally, to insure sufficiency of t.ie above conditions, recall that 
the solution of an n order differential equation generates exactly n 
linearly independent solutions. The solution for the n modes of the 
system, then, form sets of linearly independent vectors, provided all 
modes are present. 

T 

From the fact that R^ has the decomposition X X, it follows that 

it is positive semidefinite. If the system is identifiable, R will be 

o 

nonsingular and therefore positive definite. R is further seen to be 

T ° 

symmetric, from the decomposition X X. 

6 . IDENTIFICATION OF A SYSTEM WITH CONSTRAINED OBSERVABILITY 
The material presented in this section is an extension of a de> 
veloixnent of Lee [33* who considers the case of a scalar observable. 
Consider first a free (unforced) linear process of order n. 

Let the observability of the state vector be constrained 

Zk » HXj^ (6.2) 

Here z^^, the observation. Is an (ix.I) vector, Kn, and H, the 

observability transformation is of dimension (ix n) , Such a system 

is said to be observable if the initial state vector x may be deter- 

o 

mined from the sequence of observations z , z,, . . . assuming 

o 1 

that H and 4 are known. Forming this sequence of observations into 

} 

an (rucl) vector gives 
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j f 


, (6.3) 


where z denotes that observation required to make y of dimension 
u o 

(nxl)^ For example, if n»5 and 1 *2 , 




Here z denotes the element z. . 

V ^ ■ j : 

From (6.1) and (6.2) 

■i I , t ) 





z 

• 'i 

H 

J 

"h o ” 


0 


0 



1'. i 

z. 


H 4 X 


H 4 


1 


o 


f 


• 


• 


• 

* 

i j: 

~ m 


K 

1 

o 

• 


• 


• 

' ' ) 





l 


• 


• 


• 

i ‘ > 

Z'- 


H ♦f'x 


H « *' 


_ V '_ 


„ 2. 


__ P _ 


X ■AX 

o o 


J-/ 


m 


> 




' I 



, '■) 




(6.4) 


tj 


defined in the eame context as z^above>fand may be written , 
out for the example of m*5> i*2. r ' ^ .i ^ 

■i c- 


o.ta' 


.0 


H 


■ rr- j*:, 

5.?. 

B 


» 'Ij lol .,^y.or‘4 -■ j-a. H 

-- • Iv4;- . -t 







Referring to (6.4), it Is seen that the unique determination of the 

vector X from the sequence of observations z , z., z , is de- 

o o I V 

pendent upon the nonsingularity of the matrix A. If A is nonsingular, 
H, 4* are said to constitute an observable pair. 

Assume H, 4 are an observable pair, and consider the forced 


system 




Define a new vector 

where the transformation A is defined by (6.4). 
Incrementing the index gives 


(6.5) 


( 6 . 6 ) 


(6.7) 


Since A is by assumption nonsingular 

\+i‘***'*>'k**^ "k 


( 6 . 8 ) 


Define 4 * « A 4 A* 


(6.9) 


From the definition of (6.9), 4 and 4* are similar matrices, 
and will therefore have identical eigenvalues. Further, if the process 
is constrained with an observability matrix H 


H Xj^ * H A" y^ * H* y^ 


( 6 . 10 ) 


From (6.10), note the sequence of observations z , z, , . . . 
will be identical if generated by either of the processes 


''k+i***yk**^ *k 


‘k -“‘''k 


( 6 . 11 ) 


Vl'* *k*^ «k 


* H X. 
k k 


( 6 . 12 ) 
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4 * and H* can be shown to be of special form. To demonstrate, 
consider again a 5*^" order system with two system modes observable. 
That is, a (2x1) vector, and Xj^ a (5x1) vector. From (6.9), 


H yi 

H« (6.13) 


where, as before, 

, "h ♦>' 

[-V- 

Taking 4 into the left bracket of (6.13) gives 
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From equation (6.15) 


CH_j H] 


1 0 
0 1 


[C_J c] 


1 0 0 
0 1 0 
0 0 1 


• (:i ' 


r. 


[H_i C] 


0 0 
0 0 
0 0 


[C_J H] 


0 0 0 
0 0 0 


performing the multiplication in (6.14) gives 

0 0 I “ 


♦ * 


0 0 I I 

0 0 I 


® H-i! * c.j 


(6.16) 


In general, by straightforward extension of the above procedure 


I , (■.■• 


^ 4r o 




i 


‘ 0 . . .1 


T 


T 

■ 


(6.17) 


H* can be found using (6.10) and (6.15), for the example case 
of H a (2x5) matrix. ' ^ ‘ 


H* ^ H A"^ • H t <5 C_j] » 


1 ii >* ' "■ ' w ;'; i 


: f < u 


I • j 


1 0 0 0 0 

0 . 1 , 0 0 0 
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For the general case 


H* 


0 . . . 


I 


(6.18) 


where H* da of dimension (i xn). 

Reconsider equation (6.11) for the special case where the first 
(n-ji) elements of the vector A F w* are identically aero for all w*. 

•k -“*>'k 

A moment's reflection of the special form of 4 * and H* will show 

a. 


k+1 


k+v 


(6.19) 


Thus, for this special case (the conditions being on A F), the com- ^ 
plete state vector of an equivalent linear system ^ H*, may be 

constructed from the observations , provided only that H, ^ form an 
observable pair. The case of nonaero elements in the vector A Iw* 
is disdussed in Section 8. 

7. SOME PROPERTIES OF THE ESTIMATOR , 

This section is concerned with some properties of the estimator, 

^ . Recall (3.6) and (3.9), 




(3.6) 


A T T -1 
- Vj X (X^X) 


T -1 
V, X R 
i o 




(3 .'9) 




\ r. f i 

This fcxmulation fits nicely into the framework of classical re- 
gression analysis [4 ] . ^ This foot can be utiliaed to determine several 


2d 











A 

statistical properties of the estimator, , In particular, the ex¬ 
pected value of the estimator may be calculated. Recall the vector 
was formed from a scalar sequence of zero mean. 




E(w|^ = 0 


From (3.6) and (3.9) 


E(v^) « E(X + c^) * X 


(7.1) 


' E( E(vJ X (X‘X)'b - (7.2) 

If has non zero mean, say E( ) « d then 

E(ij) -0^ + 6 jXR“^ (7.3) 

A 

A requirement, then, that be an unbiased estimator is that the 
excitation be of zero mean. 

To calculate the covariance of the estimator 

.. . j . n i , 

cov( 0j) - E ( 0^ - E 0j)( 0J - E 0^)^ (7.4) 

note that ’ ■ ‘ . 

0j - E 0^ - (x'^X)"^x'’^(Vj - E v^) (7.5) 


V 


i 



cov(Vj) ■ E( I { a scalar) (7.6) 

'^1.4. j’-* - e '<1 1 'i - .1- f. 

If is formed from the sampling of a band limited spectrum, 
(7.6) will not hold, since it assumes 


V 

E(wJ^wJ) - 



J-k 
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.X’ 


2 

is further seen to be equal to q^^, where Q * cov(W|^). 

Assuming be from a white spectrum gives, using (7.5) and (7.6) 

cov( 0.) - (X R"^)^ E( C. cj) X R“^ - R"^ (7.7) 

1 O i 1 O 1 o 

For the case of a gausslan sample of zero mean, one can write 
out the Joint distribution (likelihood) fimction of the sample. Desig¬ 
nate the log of this function by L( 0^ , a^). 

L( 0j » j) = -im(log 2 + log a*) 

- i(Vj - X 0^)'^(vj - X 0^) (7.8) 

The maximum likelihood estimator is formed by setting 

Note however, from equation (3.7), that the gradient of the original 

2 

cost function g., is precisely the gradient of L( 0. , o ). Clearly, 

A * i i 

then, 0^ is the maximum likelihood estimator in the case of gaussian, 

11 

zero mean excitation. 

To find the maximum likelihood estimator for the variance of , 

2 

the excitation, a, , set 
i 1 

Replacing 0, with 0, in (7.8), and using (3.9), this expression 

1 1 1 :i .> t.r) '■ 

reduces to 

* » ■) , ' i . 

if * - 0j X)^Vj/m . (7.9) 

/ -I 
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This same estimator can be derived from the sample variance of 
since 

A "‘■1 T T 

■ k-O "k-0 

, A 

Replacing 4 with 4 , and using the fact that 

gives 

Q • S i^(w wM « (R - 4 R^ )/m (7.11) 

k*»0 ° 

This estimator must converge to Q, since 4 Is unbiased. Note 
that the same result for (q^^) Is given by (7.9). For the case of F an 
n vector and w* a scalar, (7.11) can be used to estimate F , for In 
this case, 

Q «E(Fw* wj^F^) - a ^ FF*^ (7.12) 

The properties of the estimator 4 are summarized as follows. 

A 

1) 4 is the minimum variance, unbiased estimator of 4 , provided 
the mean of the excitation is zero. 

2) For the case of gaussian, zero mean excitation, 4 is the 
IV ) maximum likelihood estimator of 4 . 

3) For the case of gaussian, zero mean excitation, 

• V • :1 -t 

cov ( 0j) -Oj R^^ 

Q « (R^ - < Rj )/m , 

8 . IDENTIFICATION WITH NOISY AND CONSTRAINED OBSERVATIONS 
Consider first a system whose complete state vector may be ob¬ 
served, but the observations are contaminated by additive noise. Let 
the observation noise be from a gaussian, zero mean distribution. 
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= ♦ X, + w, 

k k 


( 8 . 1 ) 




In consonance with the definitions of soctlbh 3, define 
i 


,, i 


1 


i 


n 


K h-- Vi> 


(nxm matrix) 


(mxl vector) 


Defining X, and as in section 3 gives 


Vj “ X 0^ + 


(3.6) 


Introducing the new variables 


“i i 




: V. 


W -X + R 


gives« from equation (3.6) 


Uj - Oj » (W - R) + Cj 


. •• . f 


(8.2) 


.£ ' ! 


Note here that the arrays u^ and W may be formed from the sequence 


of noisy observations« z^. 


Rearranging (8.2) gives 

Uj * W 0^ + - R 0^ 


^!Now define a nlew vector 




-du otJ lo : . I 

Ci“ri+Pi“R0i 

’ 9 -‘ 


t • A 




A typicel element of'this vector is given by 


\c. 


>k i ^ 1 T ^ 
«i ■ *k*'k-'k-i*i 
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From the assumed distribution of the random vectors r and w 


E(r,r*)-0 

E(w^ rj) - 0 


for J / k 
for all J, k 


Thus is the sum of independent, zero mean gaussian vectors, 
and is therefore gaussian with zero mean. The variance of the vector 
is the sum of the variances of its component parts, 

j,-;] ' 


and from (8.2) 


‘i 




(8.3) 


(8.3), and the conditions on 4 ^ it 


From the form of equation 
is apparent that all the preceeding analysis leading to an optimum iden¬ 
tification scheme is also valid for the case of noisy observations. The 
deleterious effect of the measurement noise is readily apparent in the 
increased covariance of the estimator. 

cov(0 )-a? (w'^W)"^ 

. The formulation also applies to estimation with constrained ob- 

\ • 

servations, as discussed in section 6. 


^k+1 


♦ * Yi. + A r w. 
k k 


H*y, 


Here z^^ is an ( i xl) vector, t<n. Making the assumption that 
:irst (n-i 
it was shown 


the first (n>i) elements of the vector A F w* were Identically zero. 


^k- 


{ .-) 


k+1 


z 


k+p 


(6.19) 
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For the more general case where A F w* has non zero elements, 

th ^ 

consider again the example of a 5 " order system with 2 observables. 

Let r be a vector distribution matrix, and w* a scalar. Let the 

1 2 5 

elements of the vector A F be designated d , d , . . . d . From 
equations (6.11) and the designation of A F , and utilizing the form of 
H*, , 






(8.7) 


The f pectral density of r* is clearly a function of A F . Comparing 
(8.6) and (8.1) 

( 8 . 6 ) 

*k+i ■ **k * "k 

( 8 . 1 ) 

it is noted that the format of the two sets of equations is identical. 
Thus, the net effect of forming “pseudo" state vectors from the se¬ 
quential observations of dimensionality l<n is seen, from (8.6), to 
be equivalent to adding correlated (colored) measurement noise r* to 
the state vector y^^ of an equivalent linear system, 4> * . 

The correlation products of r* produce an apparent impasse in the 
application of the identification methods herein presented. The prop- 

A 

erties of 0^ as developed in section 7 are dependent upon the assumption 
that the perturbing sequence have the property of statistical Indepen- 

A 

dence. The unbiased and minimum variance properties of then, do 
not necessarily hold in the case of constrained observability. Com¬ 
putational investigation has verified that, in general, ^ * will not be 
an unbiased estimator. 

Arriving at this same Juncture, but through a different approach, 

Lee C 3 3 suggested using only every n^ observation in forming the 
estimator 4 * . However, frcsn (8.5), it is apparent that no amount of 
time separation in the observations utilized will render r* an ta- 
correlated sequence. C<xnputatiQnal experimentation has verified this 
assertion, and is demonstrated in section 10. 

Another approach considered was to extend the order of the es¬ 
timator to include a coloring filter. Since the measurement noise is 


^k+1 




+ A r w* 
k 


''k * 'k 
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lumped with the excitation to form 4^ in (8.3), dnd since there must 
exist a linear filter or order (n-1) that will produce r* from a white 

A 

spectrum, it was reasoned that an estimator 4** of order 2n>l might 
have the necessary combined characteristics of the system being iden¬ 
tified and the coloring filter. This line of reasoning also appears to be 
invalid, however, as is demonstrated in section 10. 

9. A COMPARISON WITH OTHER ESTIMATION SCHEMES 

Ho and Walen [ 2 ] have presented the recursion fonnula 

i “XW - <*•!> 


where X is a scalar constant. Using this formulation, the con¬ 


vergence 


11m 



4 


'I' 


can be demonstrated with probability 1, provided an appropriate value 
of X Is selected. 

j 

Compare this formulation with (4.S) 




X ) x^ P 
m m m 


(4,5)" 




a *,(x P x*'^ + 1)”^ 
m m m m 


I • ■ j - 

In (9.1) one can consider the term 


(x ^ X ) x^ as being the 
' m+1 m m' m 


current estimate (based on the transition from x to x_^,) of the 

m m+i 

error in the estimator> 4 . Note that this term also appears in (4.5). 

m 

Extending this line of reasoning, the factor (1/Xm) in (9.1) is the* 

A 

weighting given this error term in-forming the new estimator, 4^^^^. 
Setting X ** 1 will^bause every error so generated to be weighted .< 
equally in the aggregate estimate. For X< 1', later terms will be 
weighted rriore heavily than the earlier terms, "nr' i:* 

•t k:*!3n'r.*i-.. li ‘‘ j u l-:r, 


28 


The 


The analogous weighting factor in (4.5) is the term a P . 

m m 

two formulations are seen to be identical except for this factor. Thus. 
(9.1) may be considered to be a special case of (4.5). with non¬ 
optimum weighting of the error term (optimality here being taken in the 
usual least squares sense). A comparison of the schemes with 1 
is presented in section 10. 

Lee [ 3 ] considers the identification of a system with a scalar 
obseivable. '’ The input distribution matrix and state transition matrix 
are constrained to be of the form 



This system is described in section 6. Lee's formulation for the iden¬ 
tification df this system is equivalent to that herein presented. 

I 

10. computational RESULTS 

The purpqse of this section is to obtain computational substan¬ 
tiation dftiiS formulations herein presented. To this end. the follow¬ 
ing experimental objectives have been set forth; 

1; '^est die algorithm 4 • R, R”^ for convergence to ♦ . 

1 o 

2.' Test the recursive algorithm for convergence. 

3^‘ Compare 1. and 2. above. 

4. ifest the recursive algorithm with additive measurement noise. 

5'.'^‘Investigate convergence under varying observability con¬ 
straints . 

‘6.^ Compare Ho'8 [ 2 ] method (equation 9.1) with the recursive 
algorithm 

Two model plants were selected to implement the testing. The 
first is a simple 'oscillator; the second a fourth order plant with two 
oscillatory modes, taken to be representative of the longitudinal 
dynamics of a large jet transport aircraft during ncxmal cruise .Cl] 





Parameters for the two system are as follows: 

2 ”^ order plant: w ■ 1 rad/sec T ■ ir/8 sec 

4^^ order plant: ■ 1.15 rad/sec ■ .35 

<02 ® .11 rad/sec ■ .035 

T «= 1.5 sec 

Testing for convergence , (tests 1, 2 and 3) 

A 

All the experimentation done supported the contention that 4> does 
in fact converge to 4 , and that the estimator ^ is therefore unbiased. 
The results of testing with the 2*^^ order oscillator are presented as 
being typical. In fig. (3), the magnitude of the elements of ^ are 
plotted as a function of the number of observations used in forming the 
estimator. These calculations were made using the nonrecursive (batch 
processing) algorithm. 

A 

In fig. (4) , the elements of 4 are plotted as a function of the 
number of iterations« and were calculated using the recursive algor- 

g 

ithm. The recursive algorithm was Initialized by setting P. * 10 • I. 

A 2 ^ 

Since it was demonstrated in section 7 that cov(0.) * a. P_. > the 

1 1 m 

large initial value for P^ demonstrates this uncertainty in the initial 
value of ^ ^, taken to be the identity matrix. The resultant estimator 
is seen to be very nearly equivalent to that of the batch processing 
technique, by a comparison of figs. (3) and (4). Using this initial¬ 
ization scheme for the example shown, all elements of the recursive 
estimator were to within 4 significant figures of the nonrecursiye 
estimator at 300 iterations. Identical data was used in the generation 
of the two estimators. 

elt is noted that the recursive form of the algorithm offers several 
computational advantages, including its suitability for real time im- 
plemeptstion and the fact that no matrix inversion is required in its 
implementation. Because of these advantages, the recursive form 
was used for the majority of the remaining testing, ^the author being 

A ' 

satisfied that no degradation of the estimator 4 would result from so 
doing. 
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Iterations 


Figure 4. 
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Estimator 0 , using the recursive algorithm, 2 


nd 


200 


order plant. 
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I 

.1 


Testing with measurement nol«e. (test 4) 


For the remaining tests, especially with the higher order system, 

It Is desirable to Introduce a scalar measure of merit for the estimators. 
While the definition of such a measure Is certainly arbitrary. It Is the 
author's opinion that such a measure should make a comparison of the 
characteristic roots of the estimator with the characteristic roots of 
the actual system. This method seems especially appropriate when 
comparing 4 * matrices, the make up of which vary considerably with 
varying observability constraints. By definition, then. 



where the X are the characteristic roots of the matrix 4 , and the 

A * A 

X ^ the roots of the estimator 4 . 

A typical result of estimating with uncorrelated measurement 
noise, using the 4^ order model. Is shown In Fig. (5). Here for 
the case of noisy observations Is compared with an estimator formed 
from the same data, but without the additive measurement noise. For 
this example, the ratio of variances for excitation and measurement 
noise was taken to be 


10 forl,J - 1, 2, 3, 4. 

Wi r^ 

Testing with constrained observability, (test 5.) 

A 

As stated In section 8, the estimators 4* , formed from systems 
having constrained observability. In general are not unbiased. How¬ 
ever, several Interesting properties of the estimator have come Into 
evidence from computational experimentation. 

The first experiment under constrained observability was con- 
ducted with the 2 order oscillator, and the estimator 4 * formed 
from observations-on a Single mode of the system. 












Iterations 


V 







0 


Figure 6. Estimator 


A 

0 * 


Iterations —^ 

V 

for a simple oscillator. Estimator formed from 


ro:> TO observations on a tingle mode of the system. 

. er< . 7af'to • , no* r -r.^au oju; 
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Results of the experiment are shown in figs. (6) and (7). Convergence 
to 4 * , or to a matrix very close to ^ , is demonstrated in this case. 

Fig. (7) shows a comparison of the constrained observability estimator. 

A A 

. and the estimator formed from the full state vector. 4 . The 
measure of merit suggested by Lee [ 3 ] 


II* -^11^ 

^' II * II ^ 

-6 

is plotted here. At 300 interations. M. equalled 10 for the case of 

-4 ^ 

the vector observable; 5x10 for the scalar observable. 

Experimentation with the 4^^ order system under constrained ob¬ 
servability failed to produce a convergent estimator. However, five 
distinct implementations had the commoniproperty of correctly iden¬ 
tifying the dominant poles of the system. Fig. (8) shows the charac¬ 
teristic roots, at 300 iterations, of 4 estimators, based respectively 
on a scalar. 2 vector. 3 vector, and full state vector observable. At 
2000 iterations, the roots of the estimator based on the full state 
vectcx had moved to within 3 significant figures of the actual roots of 


the system. while the roots of the other three estimators had not moved 
appreciably frcmi the positions plotted in fig. (8). 


Fig. (9) compares the roots of an estimator formed by using only 
th 

every 4 observation with an estimator formed|a scalar observable at 

th 

each sampling instant. Using only every 4 sample Is the method 


suggested by Lee [ 3 ] . as discussed in section 8. Note here that 

the state transition for 4 sampling instants is given] 4* (4T) « 

4 4 ^ 

[ 4 J (T) ] . The roots of 4 J (4T) are then ( X j) , where X are the 

roots of 4 f (T). Making this calculation places the dominant poles of 

4i(T} in close proximity of the poles of 4, as shown in fig. (10). 

In section 8 it was suggested that the order of the estimator 4 ** 

■ . -i ' j . . 

for the case of constrained observability be increased to 2n-l to allow 
for the inclusion of a coloring filter. A typical result of such an 


I 

augmented estimator is shown in fig. (11). Also shown, for comparison. 




Z PLANE PLOT OF CHARACTERISTIC ROOTS 



+ Roots of the actual state transition matrix 
o Roots of 4* , based on a scalar observable 
^ Roots of ♦ J , based on a 2 vector observable 
D Roots of ^ i , based on'a 3 vector observable 

A ^ \ > 

L Roots of 4 . All state modes observable 

Figure 8. The characteristic roots of 4 estimators are compared. Only 
the estimator 4 Is converging to ♦. The dominant (low fre- 

f. * 

quency) roots of all estimators are within 3 significant 
figures of the dominant roots of ♦ . 300 iterations were 

i ’■ n. ■ ■ ■ i:. f ... ■■ 

used to generate all estimators. 

•; i 'U , .'t. • -..1 -, * ■ 
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Z PLANE PLOT OF CHARACTERISTIC ROOTS 



Roots of the actual state txansttlon matrix 

O Roots of 4 , based on a scalar observable 

^ Roots of 4 i , based on a scalar observable, 

* th 

using only every 4 observation 

Figure 9, The characteristic roots of an estimator ^ S , formed from 

th ^ * 

jevery 4 observation, are compared with the roots of 4 

farmed from successive scalar observations. 2000 iterations 

were used to generate the estimators. The dominant (low 

frequency) roots of ^ ^ are within 3 significant figures of 

the roots of 4 . 




Z PLANE PLOT OF CHARACTERISTIC ROOTS 



) 

Roots of the actual state transition matrix 

A i ' 

O Roots of ♦ * 

.1 . 

^ Roots of ♦ *(T) 

« i 


Figure 101 The characteristic roots of 4*(T) ar6 generated from 
^|{4T). these roots are compared with the rdots of 
The data of fig. (9) Was u^ed to calculate these 
roots. 



ar« the rooti of a 4^ order estimator based on a scalar observable. 

2000 iterations were used to generate the data shown. 

It is of interest to note the proximity of the roots of the 4^ order 
estimator to 4 of the roots of the augmented estimator. 

Comparison with Ho's method, (test 6.) 

Figure (12) shows the location of the characteristic roots of an 
estimator generated from (9.1), the formulation proposed by Ho. [2 ] 
These roots are compared with the roots of an estimator formed from 
identical data, using (4.5), the recursive algorithm. 300 iterations 
were used to generate these estimators. In this simulation, the roots 
of the estimator generated by (9.1) had not acquired a "convergent 
tendency" at 300 iterations, in that they still were moving quite marked¬ 
ly from iteration to Iteration. In contrast, the roots of the estimator 
generated by (4.5) were confined to a small region of the z plane after 
about 50 iterations. Thus, it is not intended toi'^ply that the estimator 
generated from Ho's scheme is not convergent, but simply to compare 
the two estimators at what is considered a typical point in time. Un¬ 
fortunately, the author did not have sufficient time to accomplish a 
more rigorous investigation of (9.1). 
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Z PLANE PLOT OF CHARACTERISTIC ROOTS 


I 



Roots of the actual state transition matrix 

A 

O Roots of 4 * based on a scalar observable 

A ^ 

^ Roots of 4 ** , augmented to order 7 


Figure 11. The characteristic roots of an augmented estimator 

A 

4 **, fonned from a scalar observable, are compared 
with the roots of 4 *. 2000 iterations were used to 
generate the estimators. 
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Z PLANE PLOT OF CHARACTERISTIC ROOTS 



't* Roots of the actiial state transition matrix 
o Roots of 4 , Full state vector observable 
^ Roots of Ho's estimator 

Figure 12. The characteristic roots of an estimator proposed by Ho 

A 

are compared with the estimator 4 , formed from the full 
state vector. 300 Iterations were used to generate the 
estimators. 
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It. AStTSACT A 

j A ^east squares estimator 4 Is derived for the state transition matrix 4 

, < 

of a linear, stationary sampled data system operating In a stochastic environment. 

; • A I 

The estimator 4 is shown to be unbiased and minimum variance under the con¬ 
dition of full observability of the state vector of the system. The estimator is 

< 

also shown ito be the Maximum Likelihood Estimator for the case of the stochastic 

i i 1 

environment having Gaussian statistics. The estimation scheme Is compared 

I , * 

with two other recently published estimation schemes, both of which are shown 
to be special cases of the scheme herein inresented. 
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